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Abstract 


Mathematical models that predict performance can aid in the understanding and development of solid oxide fuel cells (SOFCs). Of course, various 
modeling approaches exist involving different length scales. In particular, very significant advances are now taking place using microscopic models 
to understand the complex composite structures of electrodes and three-phase boundaries. Ultimately these advances should lead to predictions of 
cell behavior, which at present are measured empirically and inserted into macroscopic cell models. 

In order to achieve this ambitious goal, simulation tools based on these macroscopic models must be redesigned by matching them to the complex 
microscopic phenomena, which take place at the pore scale level. As a matter of fact, the macroscopic continuum approach essentially consists 
of applying some type of homogenization technique, which properly averages the underlying microscopic phenomena for producing measurable 
quantities. Unfortunately, these quantities in the porous electrodes of fuel cells are sometimes measurable only in principle. For this reason, this 
type of approach introduces additional uncertainties into the macroscopic models, which can significantly affect the numerical results, particularly 
their generality. 

This paper is part of an ongoing effort to address the problem by following an alternative approach. The key idea is to numerically simulate the 
underlying microscopic phenomena in an effort to bring the mathematical description nearer to actual reality. In particular, some recently developed 
mesoscopic tools appear to be very promising since the microscopic approach is, in this particular case, partially included in the numerical method 
itself. In particular, the models based on the lattice Boltzmann method (LBM) treat the problem by reproducing the collisions among particles of 
the same type, among particles belonging to different species, and finally among the species and the solid obstructions. 

Recently, a model developed by the authors was proposed which, based on LBM, models the fluid flow of reactive mixtures in randomly generated 
porous media by simulating the actual coupling interaction among the species. A parallel three-dimensional numerical code was developed in order 
to implement this model and to simulate the actual microscopic structures of SOFC porous electrodes. 

In this paper, a thin anode (50 pm) of Ni-metal/YSZ-electrolyte cermet for a high-temperature electrolyte supported SOFC was considered in 
the numerical simulations. The three-dimensional anode structure was derived by a regression analysis based on the granulometry law applied to 
some microscopic pictures obtained with an electron microscope. The numerical simulations show the spatial distribution of the mass fluxes for 
the reactants and the products of the electrochemical reactions. The described technique will allow one to design new improved materials and 
structures in order to statistically optimize these fluid paths. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 
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Nomenclature 


geometric area 

lattice speed 

mutual diffusivity 

specific internal energy 

continuous single particle distribution function 
Faraday’s constant 

acceleration due to an external field 
electric current density 

kinetic forcing term 

length 

single particle mass 

molecular mass 

collisional operator 

universal gas constant 

time 

temperature 

macroscopic velocity 

microscopic velocity 


SENT RO RSH TOR PREYS 


Greek letters 

discrete step 
porosity 

relaxation frequency 
kinetic viscosity 
density 

collision time 


ADU eM & 


Subscripts and superscripts 


value at I/O bottom value 
generic species 


A generic species 
B generic species 
e equilibrium 

m mixture 

0 

o 


cogeneration systems, and the possibility of being able to con- 
sider different composition of syngas as the fuel are interesting 
technical challenges. The Department of Energy (DOE) of the 
United States a few years ago initiated a set of research projects 
(SECA—Solid State Energy Conversion Alliance) with the 
purpose of increasing the power density, reducing the manu- 
facturing costs, and encouraging commercially cost-effective 
prototypes. The European Union, through the European Hydro- 
gen and Fuel Cell Technology Platform: Strategic Research 
Agenda (January 2005), has indicated that the SOFC is a priority 
choice for stationary applications. 

The typical solid oxide structure involves a thick electrolyte 
of yttria stabilized zirconia (ZrO2 + 8% Y203) for the electrolyte 
supported (ES) cells. The function of the YSZ ceramic sup- 
port is to maintain the stability of the electronically conductive 
nickel—metal particles and to provide an anode thermal expan- 
sion coefficient acceptably close to those of the electrolyte. The 
YSZ part of the cermet structure also partially serves as the 


ionic conductor. Very thin layers (50-100 um are typical val- 
ues) are coated on both sides in order to form the cathode and 
the anode elements and must sustain actual operating conditions 
of 900-1000 °C. 

A new trend during recent years is to employ thick porous 
electrodes, which are made of cermet [4—6]. This practice allows 
one to lower the operating temperature to a moderate range 
(700-850 °C) and consequently to reduce the material cost and 
the requirements for auxiliaries. With regards to the cathode 
electrode, a material addressing the technical requirements of 
this component is lanthanum manganite suitably doped with 
alkaline (calcium) and rare earth (strontium) elements in order 
to improve its electronic conductivity (Lag,7Sro.3MnO3). This 
material allows one to reduce the cathode layer thickness to 
1 mm in cathode supported (CS) cells and 0.1 mm in anode sup- 
ported (AS) ones. With regards to the anode electrode, its basis 
material is a skeleton of yttria stabilized zirconia (YSZ) around 
Ni particles. The anode layer thickness is on the order of 0.1 mm 
(CS) and 1.5mm (AS). 

In many international research centers, most of investigations 
are focused on the experimental analysis of single cells SOFC 
under typical operating conditions or of small stacks fed from 
various fuels such as methane or hydrogen [7]. The objective of 
these experiments is to understand the phenomena that are at the 
base of the operation and the stability of the cells [8], character- 
izing the electrochemical behavior of a single cell’s components 
and identifying the optimal conditions for its operation [9,10]. 
Cell or stack performance is recorded with respect to parame- 
ters such as fuel composition, the pressure of the reactants, the 
operating temperature of the cell, eventual phenomena of macro- 
scopic degradation (e.g., de-lamination or carbon deposition 
[8]), cell polarization losses related to processes of microscopic 
degradation (e.g., evolution of the microstructure of the anode 
[11-15], cathode [16,17], and electrolyte [18—20]). Sometimes 
these phenomena are investigated by resorting to morphological 
analyses of the materials before and after the phases of opera- 
tion through the use of sophisticated techniques such as scanning 
electron microscopy (SEM). 

Some of the microscopic phenomena investigated in the lit- 
erature include those for the cathode electrode, which can be 
affected by kinetic de-mixing. This effect does not seem related 
to the cell operating temperature but instead to the fact that the 
material is affected by the electric field present. This condition 
of kinetic de-mixing can lead to pore formation in the material 
which is primarily situated at the cathode/electrolyte interface, 
leading in turn to losses in electrochemical performance. The 
applied current of the electric field creates an oxygen potential 
gradient, and it is this difference which may be the driving force 
for pore formation [16]. As to the anode electrode, it can be 
affected by nickel agglomeration, which can be explained as a 
degradation of the electrochemical three phase boundary (TPB) 
due to the sintering of metal particles and a decrease in contact 
area at the electrolyte/anode interface as well as a decrease in 
the specific surface area of Ni particles [15]. 

The mass transport phenomena throughout the electrode 
materials are key to cell performance. The thickness of the 
electrodes, their porosity (in volume fraction), the pore dimen- 
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sions, the tortuosity factor, and the pore connections are the 
main parameters affecting the mass transport of chemical species 
through the electrode/electrolyte interface [21]. The mass trans- 
port of chemical species greatly influences the concentration 
overpotentials on the anode and cathode sides (in particular, 
the cell voltage is primarily affected by cathode layer thick- 
ness [22,23]) and modifies the fuel reforming rate on the anode 
electrode surface, causing temperature non-uniformities and 
mechanical failure due to thermally induced stresses [24]. 

Because of the difficulties in realizing prototypes, an exten- 
sive simulation activity can be found in the open literature 
[25-27]. Many models are based on the assumption that the 
structure of the porous medium is isotropic so that the cermet can 
be described by three structural parameters (the ratio of poros- 
ity to tortuosity, the mean value and the standard deviation of 
the pore radii). They are used to compute the phenomenological 
coefficients which are involved in the constitutive correlations 
for the diffusion and the permeation (MTPM— mean transport 
pore model) [28]. In particular, the diffusion and the permeation 
phenomena depend on conventional molecular effects and on 
kinetic effects, since some of the pore radii are comparable to 
the mean free path of flowing fluids. In particular, for an anode- 
supported SOFC where high fuel utilization is required, a small 
error in the concentration overpotential calculation, which is 
mainly affect by mass transport inside the porous anode, may 
cause a dramatic change in its design performance [29]. Some 
experiments have been performed in order to determine the 
structural parameters for cermet materials [30]. These experi- 
ments and additional theoretical reasons induce one to conclude 
that among the structural parameters, it is the ratio of porosity to 
tortuosity and the mean value of pore radii that primarily affect 
the reaction rate [30]. However, at the moment, the suitability 
of structural parameters and of the constitutive correlations to 
describe macroscopic flow rates has not been completely verified 
for the cermet involved in SOFCs. This is due to reasons: 


e Many different microscopic topologies are possible which 
share the same values for structural parameters such as the 
macroscopic averaged properties of porosity and tortuosity. 

e There are a myriad of working conditions for porous elec- 
trodes and for this reason it is very difficult to experimentally 
verify the constitutive correlations which describe the diffu- 
sion and the permeation in all practical ranges. 


Despite the success of macroscopic modeling in understand- 
ing the complex phenomena occurring during fuel cell operation 
and in contributing to improved fuel cell designs, the actual 
microstructure of the porous layers that constitute a fuel cell 
is usually not modeled. Its effects on cell operation and per- 
formance are taken into account by considering homogeneous 
layers characterized by the macroscopic averaged parameters of 
porosity and tortuosity. While this eases the modeling efforts, 
it carries two disadvantages. Firstly, if an in situ measurement 
of such macroscopic quantities is performed, the related uncer- 
tainties are bound to affect the model results, whereas if porosity 
and tortuosity are treated as fitting parameters, there is no guar- 
antee that the true values are used since such parameters can be 


used to compensate for the inaccurate modeling of other phe- 
nomena. Secondly, it has been shown [31] that different porous 
layer microstructures, characterized by the same porosity, show 
different hydraulic characteristics. In other words, no macro- 
scopic parameter can exhaustively describe what happens at 
microscopic levels. 

In order to overcome these limitations, a novel approach to 
gas flow modeling in porous media, based on the lattice Boltz- 
mann methods (LBMs), has been utilized. LBMs are efficient 
numerical tools for investigating flow in highly complex geome- 
tries, such as porous media [32-34]. Even though traditional 
Navier-Stokes solvers could be used to describe porous media 
flow, LBMs do not require pressure—velocity decoupling or the 
resolution of a large system of algebraic equations [35,36]. They 
solve a simplified Boltzmann equation for an ensemble-averaged 
distribution of moving, interacting particles on a discrete lat- 
tice. The macroscopic quantities that describe the fluid flow 
can be calculated as moments of this distribution. Since the 
motion of particles is limited to fixed paths connecting lattice 
nodes, the resolution process needs only information about the 
nearest neighbor nodes. This feature, along with the explicit 
nature of the numerical scheme, makes LBMs very suitable for 
parallelization. 

Lattice Boltzmann models seem to be very promising for 
the analysis of reactive mixtures in porous layers [37,38]. For 
this reason, a lot of work has been performed in recent years in 
order to produce reliable lattice Boltzmann models for multi- 
component fluids and, in particular, for mixtures composed of 
miscible species [39-43]. The problem is to find a proper way, 
within the framework of a simplified kinetic model, of describing 
the interactions among different particles. Once this milestone is 
achieved, the extension of the model to reactive flows is straight- 
forward [44,45] and essentially involves additional source terms 
in the species equations which result from the reaction rate. 
Unfortunately, most existing lattice Boltzmann models for mix- 
tures are based on heuristic assumptions or prescribe too many 
constraints for setting the microscopic parameters, the end result 
of which is an idealized macroscopic description. 

The ultimate goal of the present work is that of obtaining a 
complete mesoscopic model of fluid flow and reaction in three- 
dimensional fuel cell porous media. The advantage of this would 
be that only the medium microstructure would need to be mea- 
sured (for example by means of microscopic images and/or 
tomography scans) and then cell performance could be pre- 
dicted. In the present paper, the following concrete steps towards 
the achievement of this goal are discussed: 


e Firstly, the reliability of the numerical simulations strongly 
depends on the reliability of the considered microscopic 
topology used in the simulations, i.e. if the microscopic topol- 
ogy actually reproduces the physical distribution of the solid 
phases. In this paper a three-dimensional microscopic topol- 
ogy was reconstructed from microscopic images using the 
granulometry law. The intrinsic limits of this technique are 
outlined. 

e Secondly, the reactive mixture model [46] previously con- 
sidered [47] was substituted with an updated model recently 
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proposed [48,49] which is based on the Gross and Krook 
model, in order to deal with reactive gas mixtures character- 
ized by large mass ratios and a non-fixed Schmidt number 
like those involved in fuel cell modeling. 

e Thirdly, a new semi-implicit algorithm, called SILBE scheme 
[49], was implemented in the three-dimensional code for 
reducing the computational time, thus, enabling larger com- 
putational domains. In particular, some numerical results are 
reported for the 256° case (16.8 million cells). 


In the next sections, the reconstruction of microscopic topolo- 
gies by means of the granulometry law, the essential features of 
the mathematical model and finally the numerical results will be 
discussed. 


2. Reconstruction of microscopic topologies 


First of all, the microscopic topology of the material for 
the numerical simulation must be obtained. For a wide range 
of materials, this can be achieved by computer tomography 
using X-ray absorption contrast (XCT). Unfortunately, mate- 
rials with fine structures do not yield the necessary image 
quality when using XCT. In particular, non-destructive X- 
ray computed micro-tomography directly produces 3D pore 
space at resolutions of around a micron. For SOFC applica- 
tion, this resolution is not sufficient and reconstructions from 
reliable 2D techniques such as standard and back scanning elec- 
tron microscopy (SEM/BSEM), is the only viable alternative. 
New three-dimensional imaging techniques like phase contrast 
tomography [50] relying on the use of synchrotron radiation 
overcome this problem. However, conventional microscopic 
images of cross-sections are still cheaper, faster, and require 
less effort to obtain [51]. 

With these images, proper digital image processing is needed 
in order to recognize the different solid phases and to distinguish 
them from the void fraction [51]. To do this, a number of proper 
threshold values are assumed in order to distinguish among the 
constituent materials of the porous medium so as to produce a 
truly digitalized image, i.e. an image where each point belongs 
to a single specified region (ion conducting, electron conducing, 
or pore). 

After 2D component identification, a proper reconstruction 
strategy (or a proper stereological theory, using the nomenclature 
of Ref. [51]) must be assumed in order to transform the 2D 
micrograph into a 3D structure, preserving the main features 
outlined by the experimental investigation. These main features 
can be computed by means of statistics of increasing complexity, 
limited only by the computational power available. 

A model of stochastic geometry is adapted to the microstruc- 
ture by fitting the model parameters. To this end, geometric 
characteristics of the microstructure are determined and the 
parameters of the model are chosen such that characteristic prop- 
erties of the material (e.g., porosity or fiber radius distribution) 
are represented correctly. 

There are two main techniques for dealing with this statistical 
regression: 


e The granulometry law, which is based on the simplifying 
assumption that the porous medium is made of regular grain 
with a fixed, known shape [52]. 

e Multiple-point statistics, which takes into account the 
neighboring information for producing a more reliable recon- 
struction [53]. 


Multiple-point statistics is used, based on two-dimensional 
(2D) thin sections as training images to generate 3D pore space 
representations. Assuming that the medium is isotropic, a 3D 
image can be generated that preserves typical patterns of the 
void space seen in the thin sections. The use of multiple-point 
statistics predicts the long-range connectivity of the structures 
better than the granulometry law. The three principal steps in 
this algorithm are as follows: borrowing multiple-point statis- 
tics from 2D training images; pattern reproduction; and final 
image processing for improving the local quality and consis- 
tency of the reconstruction. Since the implementation by the 
authors of multiple-point statistics is currently underway, the 
simpler granulometry law is used in this paper. 

Let us now consider the anode-side of a conventional 
electrolyte-supported planar SOFC. The experimental data and 
the microscopic pictures obtained by means of back scanning 
electron microscopy are taken from [54]. Information regarding 
the three-dimensional porous structure of the Ni-metal/YSZ- 
electrolyte cermet is obtained using the microscopic image of 
the porous anode obtained by a back scanning electronic micro- 
scope. The reference microscopic image used for computing the 
granulometry laws for both Ni and YSZ is shown in Fig. 1. 
First of all, some threshold must be considered in order to 
automatically distinguish among the two solid phases and the 
free pores for the mixture fluid flow. This conversion from 
the original analogical image (characterized by all the possi- 
ble values of the gray scale) to the final digital map, which 
allows one to recognize the nature of the materials involved, 
is effected by errors dealing with how to estimate the proper 
thresholds for each material. This is one of the main sources 


Fig. 1. Micrograph of the porous anode (bright: Ni; grey: electrolyte; black: 
pores) [54]. 
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Fig. 2. Granulometry law for Ni: relative populations with regards to the char- 
acteristic grain sizes for Ni. 


of errors in producing reliable post-processing of microscopic 
images. 

Once the digital map has been created, it is possible to apply 
the granulometry law. Essentially, one must count how many 
grains, characterized by a given shape (square) and a given size, 
exist in the digital map. The procedure used to obtain the total 
number of grains of a particular grain size of a particular species 
is as follows: 


e The whole computational domain is checked for the largest 
grain size possible. 

e Wherever the grain size is obtained, the corresponding cells 
are converted to neutral cells (i.e. are no longer considered in 
the following computations). 

e The number of grains corresponding to that particular grain 
size is counted. 

e This procedure is repeated until the total count for the least 
grain size of that species is obtained in order to completely 
fill the whole surface occupied by a given species. 


The results of the previous procedure are reported in Fig. 2 
for Ni and in Fig. 3 for YSZ. The granulometry laws for both 
materials are reasonably interpolated by exponential functions, 
namely: 


N d 
— = 0.4187 exp | —0.6433 — (1) 
N; Ni do 
i YSZ d/dO 
1\2 3 45 6 7 8 9 1011121314 
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Fig. 3. Granulometry law for YSZ: relative populations with regards to the 
characteristic grain sizes for YSZ. 


N d 
(+) = 0.2762 exp (-0.109<") (2) 
N; YSZ do 


where dp = 0.404 pm and is the minimum grain size considered 
by the statistical regression. Obviously the minimum grain size 
is dictated by the resolution of the micrograph. In the present 
case, the whole anode thickness (80 um) is described by means 
of 198 pixels in the considered micrograph. 

Once the granulometry law is applied for both the electron 
conducting grain distribution and the ion conducting grain dis- 
tribution, this information is used to generate the random porous 
structure of the portion of the anode layer where the three-phase 
boundaries are located. A code was developed to create the 
porous structure. The procedure used to generate the porous 
structure is as follows: 


e The computational domain is first filled with the largest pos- 
sible electron conducting grains (number of grains calculated 
from the granulometry law). 

e The next step is to randomly fill the remaining computational 
domain with the next largest electron conducting grains. 

e This procedure is continued till the smallest electron conduct- 
ing grains are filled up in the computational domain. 

e The above three steps are carried out for the ion conducting 
grains to fill up the computational domain. 


Some simple geometric relations must be applied in order to 
ensure that the 3D reconstructed microscopic topology has the 
porosity experimentally measured. 

The size of the reconstructed domain can obviously be differ- 
ent from that of the original training image used for calculating 
the granulometry law. Usually the size of the reconstructed 
domain depends on the available computational power. It is 
worth remembering that, in order to ensure grid independent 
results, the computational mesh must be finer than the physical 
grid used for describing the porous medium. In this paper, one 
of the goals is to find the optimal ratio between the size of the 
computational (structured) mesh and the physical grid, conven- 
tionally called refinement. For this reason, a small reconstructed 
domain has been selected for considering very large computa- 
tional refinements and, consequently, more accurate estimations 
of the fluid flow. In particular a 13 jm cube was considered filled 
by a 32? physical grid. The reconstructed cube sizes can be eval- 
uated in Fig. 1 by comparing them with the full anode thickness. 
The final result of the 3D reconstruction procedure is reported 
in Fig. 4. 

The next step is to locate the tree-phase boundaries (TPBs) 
where the electrochemical reactions take place. For locating the 
TBPs, it is not enough to simply check where an electron con- 
ducting cell (ECC), an ion conducting cell (ICC), and a fluid 
path cell (FPC) are in contact with each other. In fact, in order to 
realize an electrochemical reaction, the ECC must be in contact 
with the electron sink, i.e. the metal grid collecting the electrons 
on the anode-side of the fuel cell. This means that at least one 
connection path must exist between the collecting grid and the 
ECC in question. In a similar way, one must ensure that the ICC 
is physically connected with the electrolyte bulk and the FPC 
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Fig. 4. Reconstructed 3D microscopic topology used in the numerical simula- 
tions. The distribution of the three phases is consistent with the granulometry 
laws reported in Figs. 2 and 3 and is derived by the post-processing of the 
micrograph reported in Fig. 1. 
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is connected with the gas supply/discharge network (otherwise 
the reactants and the products of the electrochemical reaction 
would not be available). The location of the active TPBs, where 
the electrochemical reactions take place (by means of modified 
local concentrations of the species, as outlined in the section 
dealing with the boundary conditions), is performed automat- 
ically by the developed numerical code. This is done in two 
steps: 


e First of all, the code searches for all the regions of the recon- 
structed porous medium where an electron conducting cell 
(ECC), an ion conducting cell (ICC), and a fluid path cell 
(FPC) are in contact with each other by means of a non-null 
surface. These regions are labeled as potential TPBs and their 
geometric characteristics (essentially size and actual shape) 
are stored for post-processing. 

e In the next step, the numerical code verifies if (at least) 
three elementary paths exist which connect the potential TPB 
with the electrical grid, electrolyte bulk, and the gas sup- 
ply/discharge network, respectively. In this case (and only in 
this case), the potential TPB is marked as an actual TPB and 
the code takes into account the effects of the electrochemical 
reactions at the gas site close to it. 


In the next section, the details of the physical model used for 
describing the fluid flow of the reactive mixture are discussed. 


3. Mathematical model for reactive mixture flow in 
porous media 


3.1. Continuous kinetic model 


Following the derivation of the Boltzmann equation for a 
simple system with a single species, the kinetic equations for a 


simple system comprised of a mixture can be derived in a similar 
way [5-57]. Let us consider a mixture composed of only two 
types of particles labeled a and b. The two Boltzmann equations 
for the binary system are 


0 

A A w.V fa +ga Vea = Oaa + Can (3) 
Of 

“a tY Vfa+ 8e° Vvfs = Opa + Opp (4) 


where f4(X, v, t) is the continuous single particle distribution 
function for the a species, v the microscopic velocity, gą the 
acceleration due to an external field for the A species, and similar 
definitions hold for the B species. The quadratic expressions QO, a 
and Qpp are the collisional terms which describe the collisions 
among particles of the same type (self-collisions), while QAB 
and Qpa are the collisional terms due to the interactions among 
different species (cross-collisions). 

Each collision term has a well-known structure similar to 
the collision operator involved in the Boltzmann equation for a 
single species fluid [31]. The time evolution of the distribution 
function for each species is affected both by collisions with par- 
ticles of the same type and with particles of different type. These 
two phenomena are the kinetic driving forces of the equilibration 
process for the whole mixture. 

A simplified kinetic model which allows one to separately 
describe both the driving forces, as they appear in the origi- 
nal Boltzmann equations, would be desirable. Essentially, the 
key idea is to substitute the previous collisional terms with sim- 
plified ones, which are selected with a BGK-like structure. In 
particular, the BGK operator involves the difference between 
the actual value of the distribution function and the equilibrium 
one. Unfortunately, there is considerably more latitude in the 
choice of a linearization procedure in the case of a mixture than 
for a pure gas [57]. In the latter case, a local Maxwellian cen- 
tered on the (uniquely defined) macroscopic velocity is usually 
the candidate for the unperturbed component of the distribution 
function. In a mixture, however, it is possible to linearize about 
a local Maxwellian which contains the barycentric velocity or, 
alternatively, one can introduce distinct species flow veloci- 
ties and linearize about local Maxwellians which contain these 
quantities. 

Even though it could result in slightly more complicated mod- 
els, the choice of separate Maxwellians seems more general and 
better suitable for dealing with different regimes consistently. 
The model obtained is due to Hamel [58—60]. A Lattice Boltz- 
mann formulation of this model was recently proposed [46] and 
used for describing reactive mixtures [47]. 

Unfortunately, the Hamel model (and all those resulting from 
a proper linearization of it) are not completely self-consistent 
[48]. In fact, if one sums over the species equations, one does 
not exactly recover the momentum equations (as one should). 
Since the errors are negligible for fluid flows characterized by 
low Reynolds numbers, this model was considered in our pre- 
liminary simulations, because it allows one to tune the Schmidt 
number of the binary mixture (ratio between the molecular dif- 
fusivity and the mixture kinematic viscosity). Unfortunately, the 
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tuning strategy, i.e. finding the proper combination of relaxation 
parameters for achieving this goal, and the correction terms for 
compensating the discretization errors required some compli- 
cated algebra, which (unacceptably) increases the computational 
demand. 

A new lattice Boltzmann model was, thus, developed [49] for 
recovering the transport coefficients of real mixtures in a simpler 
way. The continuous (pseudo-) kinetic model used as the theo- 
retical starting point is the model due to Gross and Krook [61]. 
Using o =A, B for indicating the generic species, the simplified 
kinetic equation has the general form: 


dfo 1 e 
“ap TW Vo + 85° Wo = =- Uo = foim) (5) 
t Tm 


where Tpm is the relaxation time constant for cross-collisions and 
teen) is a local Maxwellian distribution function centered on a 
characteristic velocity for the mixture. The explicit expressions 
of this local Maxwellian is 


pe, = fala 1(v—uy 
J o(m) Ores)? 2 eo 


(6) 


where po is the generic species density, m, the particle mass, 
u the macroscopic barycentric velocity, e, the internal energy, 
and D is the number of physical dimensions. The barycentric 
velocity is defined as 


X opolo 
ye D a (7) 
Na 


where x, is the mass concentration (mass fraction) for the 
generic species. Local momentum conservation implies that the 
relaxation time constant Tm for the cross-collisions must be the 
same for all species. 

Macroscopic quantities, such as the density p(x, t), the 
macroscopic specific velocity ug(x, 1), and, consequently, the 
macroscopic barycentric velocity u(x, t) can be calculated as the 
moments of the density distribution function, i.e. 


+00 
Po(X, t) = ‘| Mo fo dv (8) 
+00 
Pog (Xx, t) = J MoV fo dv (9) 


This model is truly much simpler than the Hamel model; 
but it has only one tunable parameter, which is not enough for 
recovering both the desired molecular mixture diffusivity and 
the mixture kinematic viscosity. However, this problem can be 
fixed once a proper discretization of the microscopic velocities 
(lattice) is considered. 


3.2. Lattice kinetic model 


To solve the continuous kinetic equation (Eq. (5)), the dis- 
crete ordinate method can be applied [61,62]. According to 
this method, a set of discrete microscopic velocities vj must 
be defined for which the distribution function is evaluated. The 


generic function fi (x, t) is the single particle distribution func- 
tion evaluated for velocity v; at (x, t). In the present paper, a 
three-dimensional lattice called D3Q19, which makes use of 19 
discrete velocities (Q), is considered [63,64]. The specific lat- 
tice used in the calculation is identified by the magnitude c (all 
the lattices in the same set have the same proportions among 
the elements, but they are all scaled to the physical microscopic 
velocity c called the lattice speed). 

This assumption simplifies the development of the numerical 
code; and it is essential for the physical model, because the 
discrete lattice is used for overcoming the intrinsic constraints 
of the continuous Gross and Krook model. 

Thus, the kinetic equation, which is an integro-differential 

equation, reduces to a system of differential equations. The term 
that takes into account the effect of the external force field can 
be neglected, because we are interested in concentration driven 
flows (actually, a proper forcing term is considered in the numer- 
ical implementation for improving the accuracy of the method). 
Since the reference lattice has been defined, it is possible to write 
the operative formula in vectorial form, namely: 
H +V- V fo = Amlf yn) — fo] (10) 
where V is the matrix collecting all the lattice components (V has 
dimensions Q x D, i.e. 19 x 3) and the scalar product between 
matrices must be thought of as saturating the second index (in 
fact Vf, has dimensions 19 x 3 and consequently V-Vf, is a 
column vector 19 x 1). 

According to the Gross and Krook model on a lattice, the 
collisional matrix should be diagonal, namely Am =AmI, where 
Àm = 1/Tm. However, in order to increase the number of tunable 
parameters, the collisional matrix is assumed to be 


Am = Mp'DmMp (11) 


where Mp defines a proper orthonormal basis! for the D3Q19 
lattice and Dn is the diagonal matrix, namely: 


diag(Dm) = [0, A, Ah Ah ens A, ns AAN 


ml? “ml? “ml? ^m1> ^m1> ^m2> 
I M Mm Mm M M IV ,IV \IV 
Amy Am i Am: Am Ami Ami Ams Ami Am] (12) 


collecting the generalized relaxation frequencies for self and 
cross collisions. As will be clear later on (see the section concern- 
ing the numerical implementation), Ah controls the molecular 
diffusivity, aa , and a control the mixture kinematic and bulk 
viscosity respectively, while x and ay are free parameters 
effecting the stability of the model (usually Ae = ALY = 1). 
Since only the distribution functions for discrete microscopic 
velocities are considered, an interpolation test function must be 
adopted to calculate the macroscopic quantities. In this way, 
the previous integrals (Eqs. (8) and (9)) reduce to weighted 
summations of the considered discrete functions. The interpo- 
lation test function should be as similar to the local Maxwellian 


l Details on how to derive the orthonormal basis are reported in [49]. Even 
though this paper discusses the two-dimensional case (D2Q9 lattice), the deriva- 
tion technique for the three-dimensional case is the same. The transformation 
matrix is made of (simple) numerical constants. 
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distribution function as possible in order to easily include the 
equilibrium conditions. If we consider a low Reynolds num- 
ber flow (|ug| =u <|v| =c), which essentially means u «x c, the 
equilibrium distribution function can be linearized around the 
state at rest [41]. This assumption allows one to compute the 
lattice weights for performing the calculation of the macro- 
scopic quantities. The final result is that all the macroscopic 
(both hydrodynamic and not-conserved) moments are proper 
linear combinations of the discrete distribution functions; or, 
equivalently, a linear mapping exists between the macroscopic 
moments and the discrete distribution functions. 

Once the hydrodynamic moments are computed, it is neces- 
sary to verify that they satisfy the desired macroscopic transport 
equations. For achieving this goal, diffusive scaling [65] can 
be properly applied. There are three characteristic time scales 
in this system: the time scale Tc, which properly describes 
the collision phenomenon, i.e. O(tm/Tc) = 1; the time scale Tp, 
which properly describes the particle dynamics on the lattice, 
i.e. O[(L/c)/Tp] = 1 where L is the system size and, finally, the 
time scale Ts, which properly describes the slow fluid dynamics, 
i.e. O[(L/u)/Ts] = 1. The fast fluid dynamics (acoustic waves) is 
neglected. Since a lot of collisions are needed in order to travel 
across the system, Tc/Tp = ô, where ô is a small number (expan- 
sion parameter) is the porosity. Moreover, since u/c < 1, then 
Tp/Ts =ô and consequently Tc/Ts = 6”. Once the characteristic 
time scales are defined, the basic idea is to express the previ- 
ous equation in terms of some normalized quantities in order to 
analyze the slow fluid dynamics only. 

Applying the diffusive scaling to Eq. (10) for a binary mixture 
yields: 


PA LY. Coama) =0 (13) 
dB 

a + V - (pgug) = 0 (14) 
paua — u) = — DAVY pA (15) 
pp(Ug — u) = — DBV pg (16) 
PI + V-(pu@u) = Vp + vY: p(Vu + Vu”) an 


where D4 = e4/Al,, Dg = eg/àl, v = c° /(GA1) and 


P = PA + PB, PA = @APA, PB = €BPB (18) 


The internal energies (since this is an isothermal model, the 
concept of internal energy is not rigorously defined) can 
be expressed by means of the molecular weights, namely 
ea =RTIM; and eg=RT/Mg where T is the temperature and 
R is the universal gas constant (8.31441 J mol`! K~!). Conse- 
quently, it is possible to define an equivalent molecular weight 
for the mixture in such a way that p/p = RT/M, i.e. 


1 
~ xa/Ma +xB/Mg 


a9) 


This equivalent molecular weight for the mixture depends on the 
local mass concentrations (x4 = p4/p and xg = pg/p), because it 


is not an intrinsic property of the components. By selecting the 

relaxation frequency such as 
2 

al = M p 

m  MaMpg Dp 


(20) 


where D is the mutual diffusivity of the mixture, which is a 
molecular property due to the molecular interaction potentials, 
the species velocity difference can be rewritten as 
n 
PO 
(21) 


where ya =pa/p and yg=ppz/p are the volume/molar concen- 
trations. The previous expression is fully consistent with the 
Maxwell-Stefan macroscopic model, which is actually quite 
popular for modeling fluid flow in porous media involved in 
fuel cells. This is an important result because it proves that, in 
the continuous limit (i.e. if a sufficiently large number of par- 
ticles is considered), the developed mesoscopic model recovers 
the Maxwell—Stefan macroscopic model. This means that in each 
(sufficiently large) fluid pore, the mesoscopic model is consistent 
with the macroscopic approach and consequently the averaged 
macroscopic results derived from the proposed model depend 
only on the topology of the porous medium, which is the goal 
of our investigation. Before proceeding in this direction, it is 
worth highlighting that the macroscopic equations which derive 
from the lattice kinetic model do not involve chemical reactions, 
because there is no source term in the continuity equation. The 
electrochemical reactions realize a concentration driven flow by 
modifying the species concentrations at the three-phase bound- 
aries, while there are no chemical reactions in the bulk fluid. 
For this reason, the electrochemical reactions were simply mod- 
eled by means of proper boundary conditions with regards to 
the concentration for the fluid cells close to the three-phase 
boundaries. 


Mz —-M 
Vere yayB(Mpg Ay l 
YAYB M 


ua — Ug = — 


3.3. Numerical algorithm 


First of all, in order to develop the operative formulas for the 
numerical implementation, the lattice equilibrium distribution 
functions can be explicitly expressed as 


i 


HES 3 
Hon]; = Posi |p + zV 
SI C 


9 2 2 
-u+ (vi: uy —- -5u 22 
saV- 550] (22) 
where v; €v is the ith lattice velocity, while the two sets of 
constants, namely: 


o.. 1 1 1 1 1 1 
So = [1 — $50, T850; T850; T8S0> T890» 3630 3630 76 50> 

1 I 1 1 1 1 1 L l 

36 50> 3630 3650) 3650+ 36509 T890 3650> 3650 36 50° 


ASe kso] (23) 


and sj = Sọ (Sg = 1), are functions of the ratio sg = 3eg/c?. 

It is clear from Eq. (21) that it is possible to distinguish in the 
equilibrium distribution functions, terms which do not depend 
(fe) or depend linearly (fem) on the mixture velocity. Both 
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of these terms can be expressed by means of a linear mapping of 
the actual distribution functions (the equilibrium distributions 
depend on the moments and the latter are linear combinations 
of the actual distribution functions), namely: 


fn) = Meofs (24) 
Eim) = XoMe1 X fo 25) 
oO 


Consequently, the residual quadratic terms can be expressed as 


Eem = Ben — Eim — foim) 26) 


oO oO oO oO 


Secondly, some dimensionless coordinates can be considered 
(the so-called “lattice units”). Let us introduce a dimensionless 
space ĉ = x/L, time 7? = t/Ts, and collisional operator An = 
TcAm. In this dimensionless system, the continuous formula 
given by Eq. (10) can be approximated by the following discrete 
formula on a regular spatial mesh, namely: 


fê, X) -f ê- 1, =) 
= Anlfny@, X) — fo?, X) + ky, X)] (27) 


where kp, Ñ) is a proper forcing term, which must be intro- 
duced in order to satisfy the continuity equations with the 
discrete formula up to the second order in space and the first 
order in time (see [49] for details). The previous integration rule 
derives directly from the implicit backward Euler formula. The 
implicit formulation allows one to improve the stability of the 
scheme, i.e. for any values of the relaxation frequencies, the 
macroscopic transport coefficients cannot assume unphysical 
(negative) values. 

Unfortunately, the implicit formulation is computationally 
more demanding and, for this reason, the improved stabil- 
ity is paid by longer simulations. However, in this case, it is 
possible to develop a semi-implicit formulation, i.e. it is pos- 
sible to compute implicitly all the terms which can be solved 
once and for all at the beginning of the computations (pre- 
processing) and explicitly only the quadratic terms. The reason 
is that the quadratic terms do not affect the stability performance 
much. 

The first step is to derive an operative formula for the barycen- 
tric distribution function, which is computed by the sum of the 
species distribution functions. Substituting Eqs. (24)—(26) into 
Eq. (27) and summing over all the species yields: 


NEG, K) = Bs) fo — 1, $ — V) + Bod “fim @ — 1, $) 
o oO o 
(28) 


where 

N N -1 
Bs = [I — Am(Meo — D — AmMe1] (29) 
and Bo = BsAm. 

Once that the barycentric distribution function is known, i.e. 
>>, f.(7, X), then it is possible to solve the equation for each 
species distribution function (actually it is enough to solve N — 1 
equations if N is the number of species because the residual can 


be deduced by the barycentric conditions), i.e. 


fol, K) = Ssfo(? — 1, X — V) + Soff, ê — 1, Ñ) 


+ xa, X)Sc) “fo(?, $) (30) 


where 


I 


i al 
Ss = JE- An (Moo -D4 2-a 


—1 
nM (31) 


Sg = SsAm and Sc = 2/(2 — A} )SsÂmMe1. This equation is 
very similar to Eq. (27), with the exception of an additional 
coupling term for taking into account the effects due to the trans- 
port of mixture barycentric properties (computed in the previous 
step). The last term in the matrix given by Eq. (31) must com- 
pensate the additional forcing term kp, Ñ) introduced in the 
operative formula for satisfying the continuity equation with 
greater accuracy. 

The quadratic terms are the only ones explicitly computed. 
The generic concentration xo(Î, X) can be immediately com- 
puted at the new time step by means of f,(? — 1, X — V) since 
the collisional operator conserves the density (the density is the 
only conserved moment in this case because of the diffusion 
phenomenon). 

The outlined numerical algorithm reproduces the actual 
dynamics of mass transport phenomena: the barycentric veloc- 
ity field is solved first because it is the leading term and then 
the diffusion velocities for each species follow. All the previous 
matrices Bs, Bg, Ss, Sg and Sc depend only on the dimension- 
less relaxation frequencies: hence they can be computed once 
and for all at the beginning of the calculation, being the same 
for all the computational cells. 


3.3.1. Practical issues of implementation 

First of all, some coding tricks must be considered in order 
to increase the performance of the code executed on a sin- 
gle machine. The matrix notation is useful for explaining the 
essential features of the method and to analytically deduce the 
numerical accuracy of the scheme. However, the direct imple- 
mentation of the matrix formulation in the numerical code is 
highly inefficient. For this reason, all the previously outlined 
operations are directly implemented in order to outline the 
common terms. In particular, the common sub-expression elim- 
ination (CSE) was used in the design of the numerical code. 
Even though modern compilers automatically include this fea- 
ture, an accurate design of the code allows one to take specific 
advantage from the lattice Boltzmann formulation. For exam- 
ple, performing the collision step in moment space (defined by 
a proper linear mapping applied to conventional space) reduces 
the number of required operations. Moreover, eliminating the 
moment calculation as a separate step allows one to reduce the 
memory access time. 

Dealing with large computational domains (256° = 16.8 mil- 
lion cells or more), code parallelization is mandatory. In 
particular, the code has been developed in C++ and a free com- 
munication library has been adopted (MPICH 1.3) based on MPI 
technology [66]. The reported numerical results were obtained 
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on two cluster facilities. The first one is System X at Virginia 
Tech (VT). It is essentially a computational platform made of 
1100 dual-processor Apple XServe G5 (2200 total CPUs, each 
characterized by 2.3GHz, 4GB RAM and 80GB HD), con- 
nected by Cisco Gigabit Ethernet and Mellanox switches. The 
second facility has been recently developed at “Politecnico di 
Torino” and it is made of 64 Pentium 4 nodes (each character- 
ized by 2.8 GHz, 512 MB RAM and 40GB HD). In both cases, 
the typical number of CPUs used in the simulations was 64. 
Since an elementary decomposition was used, more computa- 
tional nodes would produce a reduction in the parallelization 
efficiency (see [47]), because the decomposed computational 
domain may have a local porosity which differs from the global 
one and this results in a load unbalance. The performance speed 
of the code was measured during the present tests by the actual 
MLUPS (millions of lattice updates per second). Introducing 
Nup (measured in MLUPS) as a useful parameter for estimat- 
ing code performance independently from the number of used 
nodes, namely: 


' €D? NgNc 

NUP yo Time G2) 
where € is the porosity, D the edge dimension of the cubic com- 
putational domain, Ng the number of fluids, Nc the number of 
collisions, and Ncpy is the number of CPUs. Typical values for 
this parameter found here are close to 0.28 MLUPS for System X 
and 0.14 MLUPS for the cluster facility at “Politecnico”. These 
values are slightly slower than other numerical codes developed 
for solving reactive mixtures in porous media. The reason is due 
to the fact that the semi-implicit formulation is 25% slower than 
the conventional explicit formulation. However, the stability is 
greater and, thus, for the present application, the stability of the 
method was privileged. 

Finally, the computational domain given in Fig. 4 is refined, 
which means that the computational mesh is denser than the 
physical grid used for defining the porous medium in order to 
ensure that the numerical results are mesh independent. Usu- 
ally the fluid flow in each physical cell used to describe the 
porous medium is solved by means of 2*=8 computational 
cells for the coarser meshes up to 8° =512 cells for the finer 
meshes. The typical domain decomposition for 64 CPUs is reg- 
ular (i.e. same subdivisions along the main axis) as reported in 
Fig. 5. 


4. Numerical results 
4.1. Preliminary considerations 


In the usual macroscopic modeling, some mass transport 
models inside the porous SOFC anode must be applied to 
estimate gas concentrations at the anode-electrolyte interface. 
However, such a mass transport model needs to be sophisti- 
cated enough so that it can take into account parameters such as 
complex functions of temperature, pressure, gas concentrations, 
and the physical properties of SOFC materials like porosity, tor- 
tuosity, and the pore size of the electrode materials [67]. The 


Hydrogen 
Water 
Electrons 


Fig. 5. Brief schematic of the considered application decomposed on 64 com- 
putational processes. The species flow directions and the main I/O bottom plane 
are outlined. 


essential problem is that all these models somehow depend on 
fitting parameters. 

In general, mass transport of components inside porous media 
can be described using either the extended Fick model (EFM) 
or the dusty gas model (DGM). Both FM and the DGM are 
mass transport equations taking into account Knudsen diffusion, 
molecular diffusion, and the effect of a finite pressure gradi- 
ent. Some other authors have eliminated the effect of Knudsen 
diffusion, using only the Stefan—Maxwell model in their mass 
transport equations (SMM). 

There are some evident degrees of freedom in selecting the 
fitting parameters concerning the porous medium microstruc- 
ture in order to recover the experimental data. Moreover, the 
application of the previous model requires some additional sim- 
plifications. For example, it is quite popular to neglect the effects 
due to the finite pressure gradient. This hypothesis seems rea- 
sonable since there is no net change in the number of moles 
in the gas phase due to the electrochemical reactions. From the 
practical point of view, this means neglecting the second terms 
on the right-hand sides of Eq. (21). Unfortunately, this assump- 
tion is not acceptable. The fact that the total number of moles 
in the gas phase does not change does not necessarily imply 
that the pressure is constant. In fact, the total mixture mass 
increases due to the electrochemical reactions (lighter hydro- 
gen particles are being substituted with heavier water particles) 
and, for this reason, a net out-coming mass flux must exist under 
steady state conditions. This is understandable if one takes into 
account that a relevant mass flux is entering by means of the 
(heavy) ions moving in the electrolyte. In other words, there is a 
relevant out-coming mixture flux and consequently a proper total 
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pressure gradient must exist to ensure the driving force of this 
flow. 

The goal of the calculations reported here is to directly solve 
the fluid flow of reactive mixtures in porous anodes. First of 
all, the transport coefficients must be tuned based on the (pure) 
molecular values, which come form kinetic theory. These input 
data must be the same as those considered by the macroscopic 
approaches in order to make the comparison more reliable. 

To do this, one must deal with a very large mass particle ratio 
(My,0/My, = 9) by essentially tuning the fraction of moving 
particles with regards to those at rest. Recalling that 


RT e (33) 
ea = — = SA— 
A Ma Ag 
oR e (34) 
eB = M; = SB 3 


it is possible to tune s4 = 1.00, and consequently c = 3685.8 m/s 
and sg=0.11. The corrective factor for the mixture equations 
(barycentric flux) is slightly more complicated, namely: 


a 35 
S=SA (za + xn) (35) 
because it depends on the local concentrations (in this case, 
s? =0.20 only at the beginning of the calculations and at the I/O 
bottom plane). 

Concerning the transport coefficients, the mutual ordinary 
diffusivity can be expressed by [68]: 


C. Ma +M 
D=— (eee (36) 
POxRBLD MaMp 


while the mixture kinematic viscosity can be expressed by 


XAVA XBVB 


"= T+ Fap(Ma/Mp)@p/xa) IF Foa(Ma/ Malo) 


Consequently the Schmidt number, defined as Sc = v/D, can be 
computed. The transport properties computed by means of Eqs. 
(32) and (33) are reported for in Table 1 for a mixture with 
10% by mass of hydrogen. This concentration is assumed as 
the initial condition for the numerical simulations and the trans- 
port coefficients are kept fixed during the simulations (the actual 
concentrations are very close to this condition). 

Taking into account the expressions derived by means of the 
asymptotic analysis, it is possible to correlate the dimension- 


Table 1 

Molecular properties for binary mixtures of H:/H2O 

T(°C) v (m? s7!) D (m? s7!) Schmidt 
500 0.000118 0.000201 0.589 
600 0.000164 0.000279 0.589 
700 0.000216 0.000367 0.589 
800 0.000273 0.000463 0.588 
900 0.000334 0.000568 0.588 

1000 0.000400 0.000681 0.588 

1100 0.000471 0.000802 0.587 


less mutual diffusivity with the first dimensionless relaxation 
frequency, i.e. 


x D 1 Ma+xpM 

D == (xaMa + xgBMpB) (38) 
c“ Te 3Amn Mpg 

or equivalently 

A M M 

3I _ FA + xpMB (39) 


a 3MgD 


For the mixture kinematic viscosity, a similar procedure can be 
considered with the difference that an additional term due to the 
artificial numerical viscosity must be included (preserving the 
fact that this coefficient must always be positive). Thus: 


i L(y 4 4m (40) 
v= 542, = zat = 
c? Te a), 2 
or equivalently 
k 1 
i= 41 


Finally, a proper set of boundary conditions must be consid- 
ered. Since the computational domain is chosen to be smaller 
than the physical thickness of the anode layer, periodic geomet- 
ric conditions are considered. This means that one must imagine 
that the computational domain is repeated in all directions in 
order to create the actual physical topology. 

As to the electrochemical reactions, they do not take place 
in the bulk fluid but are instead believed to take place near the 
(active) three-phase boundary of Ni, YSZ and the gas phase. It 
is formulated globally as 


H; (gas) + O? (YSZ) —> H20 (gas) + 2e7 (Ni) 


There is quite some controversy as to the actual pathway and 
nature of the elementary steps of the previous global reaction 
(see [69-72] for a general discussion of this topic). From the 
theoretical point of view, it is possible to consider various reac- 
tion pathways, focusing on the elementary-step charge-transfer 
reactions that are involved. Recently a tentative thermody- 
namic model of the H2/H20/Ni/YSZ three-phase boundary was 
established in order to investigate the influence of operating 
conditions (in terms the gas concentrations of both reactants 
and products) on the Ni/YSZ anode kinetics [69]. Three dif- 
ferent reaction pathways (oxygen spillover, hydrogen spillover, 
interstitial hydrogen transfer), based on five different elemen- 
tary charge-transfer reactions, are the basis of this model. All 
charge-transfer reactions show a strong and highly nonlinear 
dependence of their kinetics on gas-phase hydrogen and water 
concentration. However, the behavior is distinctly different for 
the various mechanisms. According to this model [69], the full 
set of intermediate electrochemical charge-transfer reactions, 
which are discussed in the following, are grouped into three 
different sets (oxygen spillover, hydrogen spillover, interstitial 
hydrogen transfer): 


(a) At the three-phase boundary line, oxygen ions (formally 
O?-) may hop from the YSZ surface to the Ni surface in an 
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elementary charge-transfer reaction (oxygen spillover) such 
as 


Oysz*” +[]ni > Oni + []ysz +2e7 


where the subscripts denote the surface to which the species 
is attached and [ ] indicates a free surface site. Subsequent 
chemical reactions of hydrogen and oxygen species, includ- 
ing adsorption of molecular hydrogen and desorption of 
water, take place on the Ni surface. Alternatively, spillover 
of hydroxyl ions (formally OH”) from the YSZ to the Ni 
surface is possible, namely: 


OHysz +[]ni > OHni + []ysz +e" 


(b) A different possible pathway for the global reaction is the 
spillover of hydrogen from the Ni to the YSZ surface. The 
hydrogen atoms may hop to either an oxygen ion site in an 
elementary charge-transfer reaction such as 
Hyi + Oysz~ —> OHysz~ + []ni + e7 
or to a hydroxyl site (both reactions may be active in parallel 
or consecutively): 

Hyni + OHysz —> H2Oysz + [Ini +e 
(c) Finally, at the typically high SOFC operating temperatures, 


both interstitial hydrogen atoms in the bulk Ni and interstitial 
protons in the bulk YSZ are known to be present in relatively 
high concentrations with high enough diffusivities to sup- 
port the global reaction. Charge transfer may, therefore, take 
place at the two-phase boundary of Ni and YSZ. Interstitial 
hydrogen and protons are formed via surface adsorption and 
surface/bulk exchange from hydrogen on Ni and water on 
YSZ, respectively. 


For each single step elementary reaction, the forward and 
backward reaction rate constants may be theoretically defined. 
By substituting the potential difference Ag by the activation 
overpotential 7, i.e. Nact = AY — A@eg, in the equations ruling 
the kinetics of the intermediate electrochemical charge-transfer 
reactions and expressing A@eq via the Nernst relation, after some 
algebraic manipulation, the familiar Butler—Volmer relation [73] 
is recovered, namely: 


azF (1—a)zF 
J = Jo 4 exp Rr t) — AP Rr Tact 


where J is the current density, Jo the exchange current density, 
z the number of electrons, œ the symmetry factor and F is the 
Faraday’s constant (F = 96,500 C mol~!). 

As clearly outlined in the previous section, the actual micro- 
scopic dynamics of the electrochemical reactions is quite a 
complicated phenomenon, which from the theoretical point of 
view is not entirely clear in all of its details. On the other hand, 
there is a large consensus in most of the applications concerning 


the heuristic Butler—Volmer relation. As of yet, the ion and elec- 
tron dynamics is not directly solved in the present model (the 
focus of our development has been elsewhere), even though it is 
not difficult to imagine some possible extensions of the model in 
this direction, since the required additional equations are essen- 
tially those of potential flow. Hence, it is not possible in our 
model to calculate at the moment the activation overpotential 
Nact for use in the Butler—Volmer relation. However, in princi- 
ple, this is not a problem for the mass transfer model discussed 
in this paper, because from the point of view of the mass transfer, 
the electrochemical reactions are modeled by proper boundary 
conditions. An eventual problem with this approach could be 
that due to the omitted back-coupling of the mass transfer on the 
reaction dynamics. 

Nonetheless, this last, more simplified approach of proper 
boundary conditions is adopted in the following. The practical 
effects on the electrochemical reactions is to locally modify the 
species concentration and hence to produce concentration driven 
flow. It is easy to correlate the total inlet/outlet fluxes (at the 
bottom I/O plane) with the averaged operating conditions of the 
fuel cell by means of the current density J so that 


J 
(paua) = +Ma 5 (42) 
T TEEN EA (43) 
2F 
Oe ere (44) 
2F 


where (-) indicates a surface averaged quantity at the I/O bottom 
plane (see Fig. 5). From the previous expressions, it is easy to 
prove that, for a fixed current, the net out-coming mass flux 
is 8=(18 — 2)/2 times larger than the in-coming one. For this 
reason, a total pressure gradient must exist as a driving force. 

The actual species concentrations close to the (active) three- 
phase boundaries must be such as to generate the mass fluxes 
given by Eqs. (41)-(43). Instead of using an iterative proce- 
dure (more computational demanding), the following strategy 
is suggested. Suppose that all the reactive cells behave in the 
same way (i.e. have homogeneous electrochemical reactions). 
The local species densities for the generic reactive cells can be 
expressed as 


PÀ = PA — Spd (45) 
Mpg 

oR = Pp + pp + dpm) (46) 
A 


where ôpp and dpm are freely tunable parameters, while p? 
and p% are the fixed densities at the I/O bottom plane. The 
ratio Mp/My, has been introduced in order to directly affect the 
partial pressures (in fact the densities must be divided by the 
corresponding molecular weights to obtain the partial pressure). 

The numerical simulations confirm that the relation between 
the mass fluxes and the tunable parameters used for specifying 
the density in the reaction fluid cells is linear, namely: 


(paua) — x9 (pu) = kpdpp (47) 


P. Asinari et al. / Journal of Power Sources 170 (2007) 359-375 371 


<PAua> 


Pafo 


ce 


140 
120 
100 


<Ppup> 


p PaCo 


< 


Fig. 6. Contours of the mass fluxes at the main I/O bottom plane (moving orthogonal to the picture): on the left hand side, out-coming (with regards to the porous 
anode) mass fluxes for the water produced by the electrochemical reactions (the black regions are where the flux is reversed due to drag effect), while on the right-hand 
side, in-coming mass fluxes for the consumed hydrogen. The normalization constants are pọ = 1 kg m~3 and co = 3685 ms~!. 


(48) 
(49) 


(pgup) — x% (pu) = —kpdpp 
(eu) = —ky dpm 


where kp and ky are the fitting parameters which can be mea- 
sured by means of the numerical simulations. Recalling Eqs. 
(41) and (42) and substituting them into the previous equations 


yields: 
(paua)| | ko —x4ku | | dep _ J | Ma (50) 
(0BuB) —kp —x%km | | ôæm| 2F | —Mp 
and consequently 

-1 
60D J kp —xko Ma (51) 
Som | 2F | kp (x4 -—Dku| |-Mes 


which allows one to select the freely tunable parameters based on 
the local densities for the reactive cells in order to produce mass 
fluxes coherent with a given current density. The simulations 
can be performed in three steps: 


e First, some guessed values for the freely tunable parameters 
ôĉpp and 6pm, are assumed. 

e Then the fitting parameters kp and ky (for the specific porous 
medium) are computed. 

e Finally, the actual current density (homogeneously dis- 
tributed) can be simulated by means of Eq. (51). 


In order to prove the effectiveness of the discussed tech- 
nique in simulating reactive mixtures for SOFCs, the following 
example is reported. The smallest grain used in the defini- 
tion of the porous medium reported in Fig. 4 is 0.404 um. A 
refinement factor of 8 was used in the numerical simulations, 
because it ensures mesh independent results (see the follow- 
ing section). This implies that the smallest computational cell 
is 0.00505? um?. The total size of the computational domain 
is 256° and was solved on 64 CPUs. Hence, each cube of the 
split domain reported in Fig. 5 is 64°. The number of collisions 
required to produce steady state conditions was roughly 60,000 
(15h of wall clock time on the cluster). 


The reactive boundary conditions were tuned for modeling a 
current density equal to /=0.4 Acm7?, which is a typical value 
for this technology. In Fig. 6, the contours for the mass fluxes at 
the main I/O bottom plane (moving orthogonal to the picture) are 
reported: in particular on the left-hand side, out-coming (with 
regards to the porous anode) mass fluxes for the water produced 
by the electrochemical reactions are shown (the black regions 
are where the flux is reversed due to back-step vortex effects), 
while, on the right hand side, in-coming mass fluxes for the 
consumed hydrogen are shown. 

Even though the local mass fluxes are very small, the actual 
transport coefficients (i.e. mutual diffusivity and kinematic vis- 
cosity) can induce some small (laminar) vortices after the solid 
obstructions which constitute the porous medium (back-step 
vortices). In Fig. 6 this is evident for the out-coming water 
because the plotted plane (I/O bottom plane) is located soon 
after the obstructions. However the same phenomenon holds for 
the incoming hydrogen too, even though the vortices are located 
in the porous medium bulk (when the pore size is large enough 
to allow them). 

Even though these considerations are very simple, they have 
led us to realize that the fluid paths of water and hydrogen can- 
not be exactly the same. The species fluid paths underlie the 
computation of the corresponding tortuosity. Hence, while the 
porosity is a geometric (static) concept, the tortuosity depends 
on the local fluid flows and, for this reason, it can be defined as 
a (fluid-dynamic) concept, which describes how a given species 
interacts with a given porous medium. It is then evident from 
the previous example that the tortuosity for the water and for 
the hydrogen should not be exactly the same, because these two 
species follow different paths. This concept will be clarified in 
the next section. 


4.2. Optimal refinement for fluid flow 


The optimal refinement (i.e. how many computational 
cells must be considered for modeling the fluid flow in a 
physical cell used to define the porous medium) directly 
determines the maximum size of the porous medium one 
can consider in the mesoscopic simulations. For this rea- 
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Table 2 
Volume averaged mass fluxes computed by different computational meshes (for 
the normalization values see Fig. 6) 


Refinement Mesh Species A Species B 
size (<Kpaua>>/poco) (<ppup>/poco) 
X2 64 25.57 3.35 
X4 1283 25.26 4.48 
X6 1923 24.56 4.77 
X8 256° 23.95 4.93 


son, the relevance of this parameter is discussed in this 
section. 

Four different meshes were considered by using an increas- 
ing number of computational cells (64°, 128°, 1923 and 256°). 
The volume averaged mass fluxes for the whole porous medium 
(identified by the <-> notation) are reported in Table 2. The 
evolution of the numerical error is plotted in Fig. 7. It is clear 
that the scheme is only first order in space: this is due to the con- 
sidered boundary conditions at the solid wall (i.e. bounce-back 
rule), which were adopted for simplicity, even though the max- 
imum order of the scheme in the bulk is higher (second order 
according to the asymptotic analysis [49]). 

This means that in order to recover the desired accuracy 
(<3%) the maximum physical grid size is that considered in this 
example, i.e. 327, solved by means of the finest computational 
mesh, i.e. 256° (refinement 8 x). Unfortunately, as it is evident 
by considering Fig. 1, a physical grid size of 327 does not allow 
one to consider a portion of the anode which is large enough 
to be truly representative of the whole electrode. However, our 
goal here is not to accurately solve the fluid flow but to properly 
characterize it from the macroscopic point of view. 


4.3. Optimal refinement for the computation of tortuosity 
Once the fluid flow is known, the macroscopic tortuosity can 
be computed. The basic idea is to interpolate the surface aver- 


aged values obtained by means of the mesoscopic simulations 
for both the mass fluxes and the density perturbations (difference 
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Fig. 7. Numerical errors in the computation of the volume averaged mass fluxes 
by using different computational meshes. 
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Fig. 8. Surface average mass fluxes computed directly from the code and inter- 
polated by a polynomial function (5th order). 


between actual density and the bulk density in the gas channel 
flow), as reported in Figs. 8 and 9, respectively. 

There is no physical principle which forces the surface aver- 
aged quantities to be continuous. Obviously the mesoscopic (or 
scale-level) quantities must be continuous (if we neglect rar- 
efaction effects) because of the conservation laws. However the 
surface averaging does not preserve the exact continuity, since 
it is applied to parallel surfaces (moving from the gas chan- 
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Fig. 9. Surface average density perturbations computed directly from the code 
and interpolated by a polynomial function (6th order). 
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Macroscopic Tortuosity t [-] 
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Depth z/L [-] 


Fig. 10. Surface averaged tortuosity for both species moving from the gas chan- 
nel deeply in the porous medium (solid curve is the average tortuosity for the 
barycentric momentum). 


nel deeply in the porous medium) which do not share the same 
microscopic topology. The interpolation is need for ensuring 
smooth functions, which can be easily derived. 

Applying the surface averaging operator to Eq. (15) yields: 
xl —d PA / dz 


= ——— 52 
Am) (paua) — (pau) ve) 


The tortuosity can be computed by taking into account that 


ED _ Ian) 


oe 53 
T= De T/M, OA 


The numerical results are reported in Fig. 10 for the same test 
case considered in Figs. 8 and 9. First of all, the computed tor- 
tuosities reasonably well match the practical values used for 
characterizing the SOFC anodes [6]. In particular a value close 
to three is quite common for different materials. Even though a 
lot of simplifications have been considered in the reconstruction 
of the porous medium, the fact that the measured tortuosity has 
a reasonable value is encouraging. 

Secondly according to the discussion reported at the end of 
the previous section, there is no reason why both species should 
have the same tortuosity, since they realize different fluid flows. 
In fact some differences are highlighted by Fig. 10. However 
these discrepancies are not so large, because proceeding deeply 
in the porous medium pressure gradients (and consequently mass 
fluxes) became very small and the computed tortuosities are less 
reliable. 

Considering the importance of tortuosity in the characteri- 
zation of the fluid flow through porous media because many 
simplified models use this parameter (e.g., [74—-76]), some addi- 
tional comments are needed in order to understand the physical 
meaning of the numerical results and their relevance with regards 
to fuel cell modeling. Despite the fact that it is usually considered 
a geometric parameter, the tortuosity was originally introduced 
as a kinematic property, equal to the relative average length of 
the flow path of a fluid particle from one side of the porous 
medium to the other [76]. If and only if a suitable model of 
the porous medium is chosen, for example, one consisting of 


Table 3 

Volume averaged tortuosity computed by different computational meshes 
Refinement Tortuosity (KT) Variance root (y/o) 
X2 2.1707 +0.4446 

X4 2.2096 +0.4519 

X6 2.2432 +0.4634 

X8 2.2733 +0.4752 


capillaries, then the tortuosity also becomes a purely geometric 
property, because, in this case, the fluid path is uniquely defined 
[76]. Moreover in many simplified models concerning the fluid 
flow through porous media (for example, the Kozeny theory), the 
tortuosity is introduced as an undetermined factor, i.e. a fitting 
property, which is usually justified only a posteriori by the desire 
to have an additional arbitrary parameter in order to increase the 
accuracy of the theoretical predictions and to take into account 
the mismatches between the microscopic topology assumed by 
the model and the actual one of the medium considered [76]. 
These three (sometimes conflicting) meanings of the tortuosity, 
i.e. as a kinematic, geometric or fitting property, are at the root of 
some common causes of confusion. In particular, if and only if 
a suitable microscopic topology of the porous medium or some 
degree of homogeneity is assumed, then the tortuosity becomes a 
purely geometric parameter, which is obviously independent of 
the nature of the flowing species. However, if a scale of descrip- 
tion small enough to catch the actual microscopic structure is 
adopted (as is the case with the present work), then the fluid 
paths are not uniquely defined and, in particular, they are not 
independent of the transport properties of the flowing species. 

Finally the sensitivity of the tortuosity to the mesh refine- 
ment has been verified and the numerical results are reported 
in Table 3. Fortunately the tortuosity has a small dependence 
on the mesh resolution (<5%). The reason is that the tortuosity 
comes from a ratio involving both pressure gradients and mass 
fluxes (as it is clear by considering Eqs. (52) and (53)), so it 
only slightly depends on the absolute errors for each singular 
quantity. From the physical point of view, since the tortuosity 
depends on the path of the considered species flowing in the 
porous medium, Table 3 proves that even very coarse meshes 
allow one to at least estimate the path of the species with accept- 
able accuracy. This means that larger physical domains can be 
simulated by means of coarser meshes (actually 64 times larger 
than that considered in this paper) for increasing the physical 
reliability of the numerical results. 


5. Conclusions 


It is clear that at the microscopic level, the actual fluid flow 
is far from being homogeneous and it is much more compli- 
cated than what is claimed by the averaged description used at 
the macroscopic level. As long as the local optimization of the 
materials is not an issue, the macroscopic description does not 
need to get involved in so many details. However, the optimiza- 
tion of the microscopic paths for both species would obviously 
produce an increase in macroscopic performance. For example, 
the numerical simulations show that the reconstructed portion of 
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the porous medium has a high reactive core surrounded by less 
reactive portions. It is possible to imagine (at least in principle) 
locating the solid phases in such a way to induce the suction flow 
of hydrogen in the central portion and the exhaust of water from 
the neighboring regions: in this way, the two species would not 
interfere with each other. From the practical point of view, this 
optimization strategy could be based on the direct computation 
of the actual tortuosity for both species in order to compare the 
performance of different micro-structured electrodes. 

Even though the reported results are preliminary, they allow 
one to qualitatively investigate how the electrochemical reac- 
tion affects the hydrodynamics inside the porous anode layer 
at the microscopic level. In particular, this proves that the 
usual practice of separately investigating the hydrodynamics and 
electrochemistry seems somewhat artificial and far from what 
actually happens microscopically. In order to apply the discussed 
approach to real-world problems, some improvements (which 
are all currently under development) must be still considered: 


e First of all, the identification of the structure’s material con- 
stituents, which is based on a number of proper threshold 
values, introduces some arbitrariness in the initial experimen- 
tal micrographs, which affect the consequent calculations. 
This problem can be avoided by means of a direct spec- 
troscopic method, i.e. the use of the energy dispersive 
spectroscopy (EDS) feature, which allows one to directly 
obtain the elemental composition of small objects or surfaces. 

e Secondly, the statistical procedure based on the granulom- 
etry law used to catch the physical connectivity among the 
constituents of the porous medium is not refined enough. 
Consequently, the consequent reconstruction based on this 
information is not completely reliable. The use of a two-point 
statistics, which directly takes into account the connectivity 
among the material constituents, would be highly welcomed 
and is under development. 

e Finally, the ion and electron dynamics must be included in 
the numerical model, meaning that two additional equations 
for the corresponding potential flows must be solved. This 
would allow one to include the Butler—Volmer relation for the 
active TPBs, in order to model the electrochemical reactions, 
which, in this way, would take into account the activation 
overpotential and the local concentrations. 
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